#!/bin/sh
#------------------------------------------------------------------------------
# =========                 |
# \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
#  \\    /   O peration     | Website:  https://openfoam.org
#   \\  /    A nd           | Copyright (C) 2025 OpenFOAM Foundation
#    \\/     M anipulation  |
#------------------------------------------------------------------------------
# License
#     This file is part of OpenFOAM.
#
#     OpenFOAM is free software: you can redistribute it and/or modify it
#     under the terms of the GNU General Public License as published by
#     the Free Software Foundation, either version 3 of the License, or
#     (at your option) any later version.
#
#     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
#     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
#     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
#     for more details.
#
#     You should have received a copy of the GNU General Public License
#     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
#
# Script
#     createGraphs
#
# Description
#     Generates velocity profile graphs for the venturiTube tutorial case
#
#------------------------------------------------------------------------------

! command -v gnuplot >/dev/null 2>&1 && \
    echo 'gnuplot not found - skipping graph creation' && \
    exit 1

# Calculation function
calc () {
    # Avoid locale problems with output of numbers from awk
    export LC_ALL=C

    awk "BEGIN{printf(\"%.4f\n\", $1)}"
}

### Venturi geometry ###
diameter=0.1

# Coordinates of geometry sections (/diameter)
venturiGeometry () {
    cat<<EOF
0.0  0.5
1.0  0.5
2.25 0.25
2.75 0.25
5.25 0.5
6.25 0.5
EOF
}

### Graph files and positions along venturi ###
graphs=""
positions=""

# Using positional parameters for graph file names
set -- A B C D E F

# Last time directory
time=$(foamListTimes -latestTime)

xLast=0.0
while read -r line
do
    # Graph file names
    graphs="$graphs postProcessing/graph$1/$time/line.xy"
    shift

    # Sampling locations at midpoints along venturi sections
    x="$(echo "$line" | cut -d" " -f1)"
    positions="$positions $(calc "0.5*$diameter*($xLast + $x)")"
    xLast=$x
done <<EOF
$(venturiGeometry)
EOF

### Gnuplot graph plotting ###
format () {
    cat<<EOF
set term pngcairo size 800,250
set output 'velocityProfiles.png'
set xlabel 'Position [m] + $scaleFactor*Speed [m/s]'
set xrange [0.0:0.8]
set xtics scale 0.5
set yrange [-0.06:0.06]
set format y "%.3f"
set ytics scale 0.5 0.025
set ylabel 'Radial distance [m]'
EOF
}

Ufmt='notitle w l ls 1 lw 1'
geomFmt='notitle w l lw 3 linecolor rgb "black"'

# Scaling factor for velocity profiles in graphs
scaleFactor=0.2

# '$Geom' variable for gnuplot to plot venturi geometry
venturiVar () {
        cat<<EOF
\$Geom <<EOD
$(venturiGeometry)
EOD
EOF
}

gnuplot <<EOF
    $(format)
    $(venturiVar)
    D=$diameter
    graphs="$graphs"
    graph(n) = word(graphs,n)
    positions="$positions"
    pos(n)=word(positions,n)
    plot for [i=1:words(graphs)] graph(i) u (pos(i)+($scaleFactor*\$2)):1 $Ufmt, \
         "\$Geom" u (D*\$1):(D*\$2) $geomFmt, \
         "\$Geom" u (D*\$1):(-D*\$2) $geomFmt
EOF
